  
# Load packages
  library(tidyverse)
  library(broom)
  library(lmtest)
  library(sandwich)
  library(car)
  library(stargazer)

# Calculate robust confidence intervals
  se_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Std. Error"]
  p_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Pr(>|t|)"]
  
  #### Rolling Cross Section Panel Study 2013
  
# Load data
  load("ObsStudy_00_data_AUTNES_Rolling_Cross_Section_Panel_2013.RData")
  dt <- table

  ### Party rating
  colnames(dt)[seq(79, 84, 1)] <- c("rat_SPO", "rat_OVP", "rat_FPO", "rat_BZO", "rat_GRUENE", "rat_STRONACH")
  
  # Set missings NA
  dt[, 79:84][dt[, 79:84] >= 77] <- NA
  
  # Coalition evaluations/preferred coalition
  colnames(dt)[seq(232, 236, 1)] <- c("rat_coal_OVP_SPO", "rat_coal_FPO_OVP", "rat_coal_FPO_SPO", "rat_coal_SPO_OVP_GRUENE", "rat_coal_FPO_OVP_STRONACH")
  # Set missings NA
  dt[, 232:236][dt[, 232:236] >= 77] <- NA
  
  # Coalition Likelihood
  colnames(dt)[seq(237, 241, 1)] <- c("coallik_OVP_SPO", "coallik_FPO_OVP", "coallik_FPO_SPO", "coallik_SPO_OVP_GRUENE", "coallik_FPO_OVP_STRONACH")
  dt[, 237:241][dt[, 237:241] >= 77] <- NA
  
  # Reverse-coding of coalition likelihood variable, because low values stand for "very likely/likely" in the original questionnaire 
  dt$coallik_OVP_SPO = recode(dt$coallik_OVP_SPO, "1=4; 2=3; 3=2; 4=1")
  dt$coallik_FPO_OVP = recode(dt$coallik_FPO_OVP, "1=4; 2=3; 3=2; 4=1")
  dt$coallik_FPO_SPO = recode(dt$coallik_FPO_SPO, "1=4; 2=3; 3=2; 4=1")
  dt$coallik_SPO_OVP_GRUENE = recode(dt$coallik_SPO_OVP_GRUENE, "1=4; 2=3; 3=2; 4=1")
  dt$coallik_FPO_OVP_STRONACH = recode(dt$coallik_FPO_OVP_STRONACH, "1=4; 2=3; 3=2; 4=1")
  
  # Dependent Variable: vote choice
  # anticipated vote choice when probability to vote larger 4
  # Set missings NA
  dt[, 209][dt[, 209] >= 12] <- NA

  # Postal voters vote choice
  dt$w1_q47
  # Set missings NA
  dt[, 212][dt[, 212] >= 12] <- NA

  # Anticipated vote choice when probability to vote smaller 5
  # Set missings NA
  dt[, 216][dt[, 216] >= 12] <- NA

  vote <- data.frame(dt$w1_q44, dt$w1_q47, dt$w1_q54)

  # Generate party choice variable
  dt$vote <- dt$w1_q44
  dt$vote <- ifelse(is.na(dt$vote), dt$w1_q47, dt$vote)
  dt$vote <- ifelse(is.na(dt$vote), dt$w1_q54, dt$vote)

  # Assign interpretable labels to vote variable
  dt$vote[dt$vote==1] <- "SPO"
  dt$vote[dt$vote==2] <- "OVP"
  dt$vote[dt$vote==3] <- "FPO"
  dt$vote[dt$vote==6] <- "GRUENE"
  dt$vote[dt$vote==9] <- "STRONACH"
  
  # Generate party choice columns based on vote variable
  dt$ptv_SPO <- ifelse(dt$vote == "SPO", 1, 0)
  dt$ptv_OVP <- ifelse(dt$vote == "OVP", 1, 0)
  dt$ptv_FPO <- ifelse(dt$vote == "FPO", 1, 0)
  dt$ptv_GRUENE <- ifelse(dt$vote == "GRUENE", 1, 0)
  dt$ptv_STRONACH <- ifelse(dt$vote == "STRONACH", 1, 0)
  
  # Coalition Evaluation 
  Z <- dt %>% select(contains("rat_coal")) %>%
    mutate_all(as.numeric) %>%
    mutate_all(funs(scales::rescale(.,to = c(0, 1))))
  
  rat_coal_SPO <- Z %>% select(rat_coal_SPO_OVP_GRUENE, rat_coal_OVP_SPO, rat_coal_FPO_SPO) %>% as.matrix()
  rat_coal_OVP <- Z %>% select(rat_coal_OVP_SPO, rat_coal_FPO_OVP, rat_coal_SPO_OVP_GRUENE, rat_coal_FPO_OVP_STRONACH) %>% as.matrix()
  rat_coal_FPO <- Z %>% select(rat_coal_FPO_OVP, rat_coal_FPO_SPO, rat_coal_FPO_OVP_STRONACH) %>% as.matrix()
  rat_coal_GRUENE <- Z %>% select(rat_coal_SPO_OVP_GRUENE) %>% as.matrix()
  rat_coal_STRONACH <- Z %>% select(rat_coal_FPO_OVP_STRONACH) %>% as.matrix()
  
  # Coalition likelihood
  gamma <- dt %>% select(contains("coallik_")) 
  
  
  coal_lik_SPO <- gamma %>% select(coallik_OVP_SPO, coallik_FPO_SPO, coallik_SPO_OVP_GRUENE) %>%
    mutate(sum_exp = exp(coallik_OVP_SPO) + exp(coallik_FPO_SPO) + exp(coallik_SPO_OVP_GRUENE),
           coallik_OVP_SPO = exp(coallik_OVP_SPO)/sum_exp,
           coallik_FPO_SPO = exp(coallik_FPO_SPO)/sum_exp,
           coallik_SPO_OVP_GRUENE = exp(coallik_SPO_OVP_GRUENE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_OVP <- gamma %>% select(coallik_OVP_SPO, coallik_FPO_OVP, coallik_SPO_OVP_GRUENE, coallik_FPO_OVP_STRONACH) %>%
    mutate(sum_exp = exp(coallik_OVP_SPO) + exp(coallik_FPO_OVP) + exp(coallik_SPO_OVP_GRUENE) + exp(coallik_FPO_OVP_STRONACH),
           coallik_OVP_SPO = exp(coallik_OVP_SPO) / sum_exp,
           coallik_FPO_OVP = exp(coallik_FPO_OVP) / sum_exp,
           coallik_SPO_OVP_GRUENE = exp(coallik_SPO_OVP_GRUENE) / sum_exp,
           coallik_FPO_OVP_STRONACH = exp(coallik_FPO_OVP_STRONACH) / sum_exp) %>%
    select(-sum_exp)  %>%
    as.matrix()
  
  coal_lik_FPO <- gamma %>% select(coallik_FPO_OVP, coallik_FPO_SPO, coallik_FPO_OVP_STRONACH) %>%
    mutate(sum_exp = exp(coallik_FPO_OVP) + exp(coallik_FPO_SPO) + exp(coallik_FPO_OVP_STRONACH),
           coallik_FPO_OVP = exp(coallik_FPO_OVP) / sum_exp,
           coallik_FPO_SPO = exp(coallik_FPO_SPO) / sum_exp,
           coallik_FPO_OVP_STRONACH = exp(coallik_FPO_OVP_STRONACH) / sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_GRUENE <- gamma %>% select(coallik_SPO_OVP_GRUENE) %>%  # Coalition Likelihood, has to be 1 as only one coalition option available
    mutate(sum_exp = exp(coallik_SPO_OVP_GRUENE),
           coallik_SPO_OVP_GRUENE = exp(coallik_SPO_OVP_GRUENE) / sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_STRONACH <- gamma %>% select(coallik_FPO_OVP_STRONACH) %>%  # Coalition Likelihood, has to be 1 as only one coalition option available
    mutate(sum_exp = exp(coallik_FPO_OVP_STRONACH),
           coallik_FPO_OVP_STRONACH = exp(coallik_FPO_OVP_STRONACH) / sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  # Mean-variance model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  N <- nrow(dt)
  
  M_SPO <- sapply(1:N, function(i) EV(coal_lik_SPO[i,], rat_coal_SPO[i,]))
  V_SPO <- sapply(1:N, function(i) Vari(coal_lik_SPO[i,], rat_coal_SPO[i,]))
  
  M_OVP <- sapply(1:N, function(i) EV(coal_lik_OVP[i,], rat_coal_OVP[i,]))
  V_OVP <- sapply(1:N, function(i) Vari(coal_lik_OVP[i,], rat_coal_OVP[i,]))
  
  M_FPO <- sapply(1:N, function(i) EV(coal_lik_FPO[i,], rat_coal_FPO[i,]))
  V_FPO <- sapply(1:N, function(i) Vari(coal_lik_FPO[i,], rat_coal_FPO[i,]))
  
  M_GRUENE <- sapply(1:N, function(i) EV(coal_lik_GRUENE[i,], rat_coal_GRUENE[i,]))
  V_GRUENE <- sapply(1:N, function(i) Vari(coal_lik_GRUENE[i,], rat_coal_GRUENE[i,])) # is 0
  
  M_STRONACH <- sapply(1:N, function(i) EV(coal_lik_STRONACH[i,], rat_coal_STRONACH[i,]))
  V_STRONACH <- sapply(1:N, function(i) Vari(coal_lik_STRONACH[i,], rat_coal_STRONACH[i,])) # is 0
  
  
  # Gender variable, 1 = male, 2 = female
  dt$male <- ifelse(dt$w1_q4 == 2, 0, dt$w1_q4)

  # Education
  # Set missings or observations that do not fit the scheme to NA
  dt$w1_q73 <- ifelse(dt$w1_q73 > 15, NA, dt$w1_q73)

  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt$w1_q31[dt$w1_q30==2] <- 0
  dt$pid_SPO <- ifelse(dt$w1_q31==1, 1, 0)
  dt$pid_OVP <- ifelse(dt$w1_q31==2, 1, 0)
  dt$pid_FPO <- ifelse(dt$w1_q31==3, 1, 0)
  
  d <- data.frame(
    "ptv_SPO" = dt$ptv_SPO,
    "ptv_OVP" = dt$ptv_OVP,
    "ptv_FPO" = dt$ptv_FPO,
    "ptv_GRUENE" = dt$ptv_GRUENE,
    "ptv_STRONACH" = dt$ptv_STRONACH,
    "rat.SPO" = dt$rat_SPO,
    "rat.OVP" = dt$rat_OVP,
    "rat.FPO" = dt$rat_FPO,
    "rat.GRUENE" = dt$rat_GRUENE,
    "rat.STRONACH" = dt$rat_STRONACH,
    "pid_SPO" = dt$pid_SPO,
    "pid_OVP" = dt$pid_OVP,
    "pid_FPO" = dt$pid_FPO,
    "lotterymean.SPO" = M_SPO,
    "lotteryvariance.SPO" = V_SPO,
    "lotterymean.OVP" = M_OVP,
    "lotteryvariance.OVP" = V_OVP,
    "lotterymean.FPO" = M_FPO,
    "lotteryvariance.FPO" = V_FPO,
    "lotterymean.GRUENE" = M_GRUENE,
    "lotteryvariance.GRUENE" = V_GRUENE,
    "lotterymean.STRONACH" = M_STRONACH,
    "lotteryvariance.STRONACH" = V_STRONACH,
    "sex" = dt$male,
    "age" = dt$alter,
    "edu" = dt$w1_q73
  )
  
  d$lotteryvariance.SPO <- scales::rescale(d$lotteryvariance.SPO, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.OVP <- scales::rescale(d$lotteryvariance.OVP, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.FPO <- scales::rescale(d$lotteryvariance.FPO, to = c(0, 1), from = c(0,0.25))
  
# Save model results
  summary(m1 <- lm(as.formula(paste("ptv_SPO ~ lotteryvariance.SPO + lotterymean.SPO + rat.SPO + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_SPO" else "")), d))
  summary(m2 <- lm(as.formula(paste("ptv_OVP ~ lotteryvariance.OVP + lotterymean.OVP + rat.OVP + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_OVP" else "")), d))
  summary(m3 <- lm(as.formula(paste("ptv_FPO ~ lotteryvariance.FPO + lotterymean.FPO + rat.FPO + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_FPO" else "")), d))
  
  AUT_2013_AUTNES_RCSP_SPO_Variance_Estimate <- tidy(m1) %>% filter(term=="lotteryvariance.SPO") %>% select("estimate")
  AUT_2013_AUTNES_RCSP_SPO_Variance_SE <- se_robust(m1)["lotteryvariance.SPO"] #tidy(m1) %>% filter(term=="lotteryvariance.SPO") %>% select("std.error")
  AUT_2013_AUTNES_RCSP_SPO_Mean_Estimate <- tidy(m1) %>% filter(term=="lotterymean.SPO") %>% select("estimate")
  AUT_2013_AUTNES_RCSP_SPO_Mean_SE <- se_robust(m1)["lotterymean.SPO"] #tidy(m1) %>% filter(term=="lotterymean.SPO") %>% select("std.error")
  
  AUT_2013_AUTNES_RCSP_OVP_Variance_Estimate <- tidy(m2) %>% filter(term=="lotteryvariance.OVP") %>% select("estimate")
  AUT_2013_AUTNES_RCSP_OVP_Variance_SE <- se_robust(m2)["lotteryvariance.OVP"] #tidy(m2) %>% filter(term=="lotteryvariance.OVP") %>% select("std.error")
  AUT_2013_AUTNES_RCSP_OVP_Mean_Estimate <- tidy(m2) %>% filter(term=="lotterymean.OVP") %>% select("estimate")
  AUT_2013_AUTNES_RCSP_OVP_Mean_SE <- se_robust(m2)["lotterymean.OVP"] #tidy(m2) %>% filter(term=="lotterymean.OVP") %>% select("std.error")
  
  AUT_2013_AUTNES_RCSP_FPO_Variance_Estimate <- tidy(m3) %>% filter(term=="lotteryvariance.FPO") %>% select("estimate")
  AUT_2013_AUTNES_RCSP_FPO_Variance_SE <- se_robust(m3)["lotteryvariance.FPO"] #tidy(m3) %>% filter(term=="lotteryvariance.FPO") %>% select("std.error")
  AUT_2013_AUTNES_RCSP_FPO_Mean_Estimate <- tidy(m3) %>% filter(term=="lotterymean.FPO") %>% select("estimate")
  AUT_2013_AUTNES_RCSP_FPO_Mean_SE <- se_robust(m3)["lotterymean.FPO"] #tidy(m3) %>% filter(term=="lotterymean.FPO") %>% select("std.error")
  
  # Harmonize names of the IVs of interest in the models
  names(m1$coefficients)[names(m1$coefficients) == "lotteryvariance.SPO"] <- "Government Lottery Variance"
  names(m1$coefficients)[names(m1$coefficients) == "lotterymean.SPO"] <- "Government Lottery Mean"
  names(m2$coefficients)[names(m2$coefficients) == "lotteryvariance.OVP"] <- "Government Lottery Variance"
  names(m2$coefficients)[names(m2$coefficients) == "lotterymean.OVP"] <- "Government Lottery Mean"
  names(m3$coefficients)[names(m3$coefficients) == "lotteryvariance.FPO"] <- "Government Lottery Variance"
  names(m3$coefficients)[names(m3$coefficients) == "lotterymean.FPO"] <- "Government Lottery Mean"
  
  
  if(!robustnesscheck_pid){
    stargazer(m1, m2, m3, title="Linear regressions of vote choice on perceived government lottery variance and mean (AUTNES Rolling Cross Section Panel Study 2013).", 
              align=TRUE, 
              dep.var.labels=c("SPO", "OVP", "FPO"),
              omit.stat=c("LL","ser","f"),
              se = lapply(list(m1, m2, m3), se_robust),
              p = lapply(list(m1, m2, m3), p_robust),
              out = "TableSM2.tex"
    )
    
  }
 
